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Abstract. We present a new way of modeling deflagration 
fronts in reactive fluids, the main emphasis being on turbulent 
thermonuclear deflagration fronts in white dwarfs undergoing a 
Type la supernova explosion. Our approach is based on a level 
set method which treats the front as a mathematical discon- 
tinuity and allows full coupling between the front geometry 
and the flow field ( |Smiljanovski et al. 1997 ). With only mi- 
nor modifications, this method can also be applied to describe 
contact discontinuities. Two different implementations are de- 
scribed and their physically correct behaviour for simple test- 
cases is shown. First results of the method applied to the con- 
crete problems of Type la supernovae and chemical hydrogen 
combustion are briefly discussed; a mor e extensive analysis o f 
our astrophysical simulations is given in Reinecke et al. (1998). 
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1. Introduction 

Numerical simulations of turbulent combustion have always 
been a challenge, mainly because of the large range of length 
scales involved. In astrophysics, prominent examples are Type 
la supernovae, where the length scales of relevant physical pro- 
cesses range from 10 _4 cm to several 10 8 cm). In the currently 
favoured scenario the explosion starts as a deflagration in the 
flamelet regime near the center of the star. At the correspond- 
ing densiti es, the typical width of the conductive flame is less 
than 1mm ( Timmes & Woosley 1992 ). Ray leigh-Tay lor unsta- 
ble blobs of hot burnt material are thought to form which rise 
and lead to shear-induced turbulence at their interface with the 
unburnt gas. This turbulence increases the effective surface area 
of the flamelets and thereby the rate of fuel consumption over 
its laminar value; the hope is that finally a fast deflagration 
might result, in agreement with phenomenological models of 
Type la explosions (Nomoto et al. 1984). 

A multidimensional direct numerical simulation of such an 
event is - and will always be - computationally infeasible; 



therefore, small scale effects like turbulence, diffusion and heat 
conduction need to be incorporated in form of phenomenolog- 
ical models. Despite considerable progress in the field of mod- 
eling turbulent combustion for astrophysical flows (see, e.g., 
Niemeyer 1995 ), the correct numerical representation of the 
thermonuclear deflagration front is still a weakness of Type la 
simulations; this is mainly due to the fact that in those simula- 
tions the conductive flame is not properly resolved, but must be 
made several orders of magnitude thicker than in reality. The 
artificially increased width of the reaction z one is a prerequi - 
site for the reactive-diffusive flame model (Khokhlov 1993), 



which has been used by most authors so far. In this approach 
the burning region is stretched out over several grid zones to 
ensure an isotropic flame propagation speed. Typical values for 
the numerical flam e width range fr om 4-5 ( Khokhlov 1993 ) 
to 8-10 grid cells ( Niemeyer 1994 ). However, the artificially 
soft transition from fuel to ashes stabilizes the front against hy- 
drodynamical instabilities on small length scales, which in turn 
results in an underestimation of the flame surface area and - 
consequently - of the total energy generation rate. 

The front tracking method described in this paper is based 
on the so-called level set technique that has been in use for 
several years in the engineering sciences. It was introduced by 
psher & Sethian (1988 ) who used the zero level set of a n- 
dimensional scalar function to represent (n — 1) -dimensional 
front geometries. Sussman et al. (1994) give equations for the 



time evolution of such a level set which is passively advected 
by a flow field; this approach can be used to track contact dis- 



continuities, for example. Smiljanovski et al. (1997) extend this 
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method to allow the tracking of fronts additionally propagat- 
ing normal to themselves, e.g. deflagrations and detonations. In 
contrast to the artificial broadening of the flame in the reaction- 
diffusion-approach, their algorithm is able to treat the front as 
an exact hydrodynamical discontinuity. Considering the fact 
that the real width of the conductive flame in a Type la su- 
pernova is several orders of magnitude smaller than the typical 
grid cell sizes in multidimensional simulations, this is a very 
good approximation. 

The outline of this paper is as follows: In section || we 
present the main ideas and governing equations of our ap- 
proach. Two different implementations of the flame model are 
described in detail in section ||. Section^ is dedicated to the re- 
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suits of simple testcases, whereas section |] lists some results of 
the application of our numerical scheme to "real world" prob- 
lems. Finally, we give a summary of open issues and an outlook 
on future work in section || 

2. The level set method 

The central aspect of our front tracking method is the associ- 
ation of the front geometry (a time-dependent set of points T) 
with an isoline of a so-called level set function G: 



r := {r | G(r) =0} 



(1) 



Since G is not completely determined by this equation, we can 
additionally postulate that G be negative in the unburnt and 
positive in the burnt regions, and that G be a "smooth" func- 
tion, which is convenient from a numerical point of view. This 
smoothness can be achieved, for example, by the additional 
constraint that 



IVGI = 1 



(2) 



in the whole computational domain, with the exception of pos- 
sible extrema and kinks of G. The ensemble of these conditions 
produces a G which is a signed distance function, i.e. the abso- 
lute value of G at any point equals the minimal front distance. 
The normal vector to the front is defined as 



VG 

Wg\ 



(3) 



and thus points towards the unburnt material. 

The task is now to find an equation for the temporal evolu- 
tion of G such that the zero level set of G behaves exactly as the 
flame. Such an expression can be obtained by the consideration 
that the total velocity of the front consists of two independent 
contributions: it is advected by the fluid motions at a speed v 
and it propagates normal to itself with a burning speed s. 

Since for deflagration waves a velocity jump usually occurs 
between the pre-front and post-front states, we must explicitly 
specify which state v and s refer to; traditionally, the values 
for the unburnt state are chosen. Therefore, one obtains for the 
total front motion 



Df = v u + s u n. 



(4) 



The total temporal derivative of G at a point P attached to the 
front must vanish, since G is, by definition, always at the 
front: 

dG P 8G _„ . dG _ n ^ 

—— = — + VG-x P = — + DfVG = (5) 

dt dt dt 

This leads to the desired differential equation describing the 
time evolution of G: 



(6) 



This equation, however, cannot be applied on the whole 
computational domain: Firstly, D f has a physical meaning in 




Fig. 1. Illustration of the basic principles o f the level set 
method according to Smiljanovski et al. (1997 ): The piece wise 
linear front cuts the mixed cells into burnt and unburnt parts, a 
is the unburnt volume fraction of a cell, (3 is the unburnt area 
fraction of a cell interface. The fluxes F u and Fb are calculated 
from the reconstructed states. 

the immediate vicinity of the front only and may be undefined 
elsewhere. Secondly, using this equation everywhere will in 
most cases destroy G's distance function property (eq. ^). As a 
consequence, this might lead to the buildup of very steep slopes 
in G which are likely to cause numerical problems (Sussman 
et al. 1994). Therefore additional measures must be taken in the 
regions away from the front to ensure a "well-behaved" | VG| 
(for implementation details, see section 



3.1.2) 



The situation is further complicated by the fact that the 
quantities v u and s u which are needed to determine D f are 
not readily available in the cells cut by the front. In a finite 
volume context, these cells contain a mixture of pre- and post- 
front states instead. Nevertheless one can assume that the con- 
served quantities (mass, momentum and total energy) of the 
mixed state satisfy the following conditions: 



p = ap u + (1 - a)p b 
~pv = ap u v u + (1 - a)pbVb 
~pe = ap u e u + (1 - ot)pb&b 



(7) 
(8) 
(9) 



Here a denotes the volume fraction of the cell occupied by the 
unburnt state. In order to reconstruct the states before and be- 
hind the flame, a nonlinear system consisting of the equations 
above, the Rankine-Hugoniot jump conditions and a burning 
rate law must be solved. The technical details are described in 



section 3.2.2 



Having obtained the reconstructed pre- and post-front states 
in the mixed cells, it is not only possible to determine Df, 
but also to separately calculate the fluxes of burnt and unburnt 
material over the cell interfaces. Consequently, the total flux 
over an interface can be expressed as a linear combination of 
burnt and unburnt fluxes weighted by the unburnt interface area 
fraction /?: 



F = (3F U + (1 - (3)F b 



(10) 



(see Fig. [[]). 



3. Implementation 

In this section we concentrate on the case of a deflagration 
wave, but the modifications needed to model contact disconti- 
nuities are straightforward: for this case, the front propagation 
speed (s or s u ) and the formation enthalpy (q) have to be zero 
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in all following equations, which leads to an overall simplifica- 
tion of the numerical scheme. 

For our calculations, the front tracking algorithm was im- 
plemented as an additional module for the hydrodynamics code 
PROMETHEUS ( fryxell et al. 1989| ). Two independent and 
completely different implementations were realized: 

- In the simpler approach, the G-function plays a somehow 
passive role: It is advected by the fluid motions and by burn- 
ing and is only used to determine the source terms for the 
reactive Euler equations. We will refer to this algorithm as 
passive implementation. It must be noted that there exists 
no real discontinuity between fuel and ashes in this case; 
the transition is smeared out over about three grid cells by 
the hydrodynamical scheme, and the level set only indicates 
where the thin flame front should be. However, the numer- 
ical flame is still considerably thinner than in the reaction- 
diffusion approach. 

- The second implementation (called complete implementa- 
tion) contai ns in-cell-reconstruction a nd flux-splitting as 
proposed by Smiljanovski et al. (1997 ); therefore it should 
exactly describe the coupling between the flame and the hy- 
drodynamic flow. 

3.1. Passive Implementation 
3.1.1. G-Transport 

Since the front motion consists of two distinct contributions, it 
is appropriate to use an operator splitting approach for the time 
evolution of G. The advection term due to the fluid velocity v p 
can be written as 



dG 
dt 



— = -v F VG, 



or in conservative form 



dt 



-v FP Gdf = 



(11) 



(12) 



( Mulder et al. 1992 ). This equation is identical to the advection 
equation of a passive scalar, like the concentration of an inert 
chemical species. Consequently, this contribution to the front 
propagation can be calculated by PROMETHEUS itself with- 
out requiring complicated modifications. As a consequence, the 
discrete values of the level set function have to be stored at the 
centers of the grid cells, like the hydrodynamical variables p, 
T, etc. 

The additional flame propagation due to burning is calcu- 
lated at the end of each time step according to the following 
procedure: 

First the four discrete spatial derivatives of G are obtained 
in each cell: 



D 



G, 



G l 



D tiJ 



D 



G, 



Gi-i 



Xi+l 

Gij+i 



Xi 



Xi 



D 



G 



Xi-l 

Gi,j-i 



Vj -Vj-l 



At the boundaries of the computational domain some of the 
above equations cannot be applied (e.g. D~ x .). In these cases, 
the gradient is set to for reflecting boundaries and extrapo- 
lated in zeroeth order for outflow boundaries. 

Afterwards, the relevant derivatives are determined by sim- 
ple upwinding with respect to the propagation direction of the 
front: 



D, 



Kii for£+ y >0andD- y . >0 
D x,ij for D tij < and /), , ; < 



(15) 



for (D 



■ D + .. 

X,IJ X,IJ 



)<o 



\D1 



where D x , tj := 0.5(1^, ,..,,„ 
The new G-value is then defined by 



Gl = G, 



3.1.2. Re-Initialization 



Ats, 



D 2 .. 



D 2 ... 



(16) 



As was mentioned in section ||, an additional correction step 
has to be applied in the regions away from the front in order 
to keep G a signed dista nce function. This task can be accom- 
plished in several ways. |Sussman et al. (1994 ), for example, 
suggest a pseudo-time approach, where the equation 



dG 
8t 



G 



\G\ 



(1-|VG|) 



(17) 



is solved iteratively until convergence is obtained. Here, e de- 
notes an empirical quantity with a value comparable to the 
length of a grid cell. While being quite efficient, this method 
has the drawback that it changes all G-values, even those near 
the front; consequently, the front might be moved by small 
amounts during the re-initialization (Sethian 1996). 

This potential problem is avoided by the following algo- 
rithm which we used for our simulations: 

- The coordinates of all zero crossings of G between neigh- 
bouring grid points are calculated by linear interpolation; 
if, e.g., Gi.j > and Gj+i j < 0, the zero crossing is at 



Gi 



G 



G 



Vz 



Vj 



{x i+ i - Xi) and (18) 
(19) 



The ensemble of all points (x z , y z ) is a discrete representa- 
tion of the zero level set. 
- For all grid points, the minimum distance to one of the 
points (x z , y z ) is determined: 



i\f(x~ l 



2 + (yj - Vz^n) 2 



(20) 



(13) - The corrected value for Gy is a weighted average of the 
original value and dij , such that 

1 ' 4) G' : Uui.X;,, + (1 - llul.j.spv.C, ,)(!,,. (21) 
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H denotes a function, which is essentially 1 for small ar- 
guments and smoothly drops to near a given threshold. In 
this work, we used the expression 



+ , 



H(d) 



1 — tanh 



d — do 



l-tanh^).(22) 



For this equation, the transition takes place in a region of 
the width 5 around do- Satisfying results have been ob- 
tained for do ~ 3A and S « A, where A represents the 
width of a grid cell. 

The weighting with H (d) has the effect that values near the 
flame are practically left unchanged, while the values farther 
away represent a distance function in good approximation. 

3.1.3. Source terms 

After the update of the level set function in each time step, the 
change of chemical composition and total energy due to burn- 
ing is calculated in the cells cut by the front. In order to obtain 
these values, the volume fraction a occupied by the unburnt 
material is determined in those cells by the following approach: 
from the value Gij and the two steepest gradients of G towards 
the front in x- and y-direction a first-order approximation G of 
the level set function is calculated; then the area fraction of cell 
ij where G < can be found easily. Based on these results, the 
new concentrations of fuel, ashes and energy are obtained: 



yi 

-* 1 Ashes 
^Fuel 



max(l 

1 

e tot + q(X' Ashes 



a, ^Ashes) 



yi 

A Ashes 



s) 



(23) 
(24) 
(25) 



In principle this means that all fuel found behind the front is 
converted to ashes and the appropriate amount of energy is re- 
leased. The maximum operator in eq. (^3j) ensures that no "re- 
verse burning" (i.e. conversion from ashes to fuel) takes place 
in the cases where the average ash concentration is higher than 
the burnt volume fraction; such a situation can occur in a few 
rare cases because of unavoidable discretization errors of the 
numerical scheme. 

3.2. Complete Implementation 

In this approach the discrete values of G are defined on the 
cell corners instead of the cell centers, because this simplifies 
the calculation of the geometrical quantities a and j3, which 
are needed for the reconstruction and flux-splitting steps. In 
the following sections all quantities defined on cell corners are 
described by fractional indices: e.g. Gi+1/2,3+1/2 denotes the 
G value in the top right corner of cell ij. 

3.2.1. Geometrical quantities 

The knowledge of the front normal n and the unburnt volume 
fraction a in the mixed cells is a prerequisite for the reconstruc- 
tion of burnt and unburnt hydrodynamical states. The normal is 



+ +, 



+ +, 



- +, 




+ 

Fig. 2. Determination of a in a mixed cell. The signs on the 
cell corners denote the sign of the G function. The rightmost 
sketch shows a situation where two different geometrical inter- 
pretations are possible and a is not uniquely defined. 

derived from the discrete gradient 



dG 
dx 



dG 
dy 



G, 


+ 1/2,3 + 1/2 


-Gi. 


-1/2,3 


fl/2 




x i+l/2 




1/2 




G, 


+1/2,3-1/2 




-1/2,3 


-1/2 




3^+1/2 




1/2 




G, 


+1/2,3 + 1/2 


-G^ 


H/2,3 


-1/2 




Vj+1/2 


~Vj- 


1/2 




G, 


-1/2J+1/2 


-G t . 


-1/2,3 


-1/2 




2/i+l/2 


-Vj- 


1/2 





+ 



According to eq. (|3j), ray is then given by 

_ _ (VG)ij 



(VG) 



(26) 



(27) 



(28) 



The value for a is found by determining the zeros of G on 
all cell edges, connecting them with straight lines and calcu- 
lating the surface area behind this approximated flame. Fig. || 
shows all topologically different situations. While calculating 
a in the first three cases is trivial, the fourth case is ambigu- 
ous since two different front geometries are possible; for this 
situation, we set a to the mean value of the two possibilities. 
Fortunately, such a geometrical constellation is quite rare in 
hydrodynamical simulations. 

3.2.2. Reconstruction 

In order to obtain the hydrodynamical state vectors U u and Ub 
from the average U in the mixed cells, a nonlinear equation 
system has to be solved. The first three equations have already 
been presented in section || (eqs. [7] - ^|). It is convenient to split 
the velocity vector into a normal and a tangential part with re- 
spect to the front; eq. (0) then reads 



pv n = ap u v ntU + (1 - a)pbV n ,b 
v t ,u = V t ,b = vt- 



and 



(29) 
(30) 



Further, the reconstructed states must satisfy the Rayleigh 
criterion and the Hugoniot jump condition for the internal en- 
ergy: 



(PuSu) 2 = - 



Pb -Pu 



= q 



Vb-Vu 

(Pb+Pu) 



(V b - V u ) 



(31) 
(32) 
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Here, e, is defined as e tot — v 2 /2 and V 
sures are given by the equation of state: 

p u =PEos(/Ou,e iitt ,X„) and 
Pb =PEOs(p&,ei,6,X b ) 



l/p. The pres- 

(33) 
(34) 



Additionally, the jump condition for the normal velocity 
component reads 



V n h 



1 



Pu_ 

Pb 



(35) 



To complete the system, a burning rate law is required. 
Usually this will be the equation for the laminar burning 
speed, depending on the unburnt state variables. In our case 
of highly turbulent burning in the flamelet regime, the flame 
spe ed can be derived from the turbul ent kinetic sub-grid energy 
e sg dNiemeyer & Hillebrandt 1995a[ |b|): 



(36) 



The ensemble of all the equations above can be solved with 
any of standard iterative method. Our implementation uses a 
globally converging Broyden solver ( Broyden 1965 ; Press et al. 



1992). In contrast to the popular Newton-Raphson approach, 
this algorithm converges even for relatively bad initial guesses, 
which is important for our application. 

3.2.3. Transport 

The algorithms presented in the three following subsections are 
designed for use with a directional splitting scheme and are 
thus orientation independent. Therefore we will only describe 
the numerical procedure for the x-sweeps. 

For the complete implementation a simple, non-conserva- 
tive approach is used to obtain the G-values at the new time 
level: 



G n 



AtDl 



dG n 
dx 



(37) 



Several complications arise from the fact that values for D x , 
which is defined in the center of the mixed cells, are needed 
at the cell corners. Since D only has a physical meaning in the 
mixed cells, its value in all other cells may be chosen arbitrarily. 
It can be shown analytically that the distance function property 
of G is preserved if the condition 



nV{Dn) = 



(38) 



is satisfied, i.e. if the flame propagation velocity is constant 
along the "field lines" of G. Consequently, the values for D in 
the whole computational domain are obtained by spreading out 
the values in the mixed cells in the direction of n and — n. 

In the next step, D x in the middle of the cell interfaces is 
calculated by simple averaging 



D 



x,i,j+l 



(39) 



burnt part 




unburnt part 



Fig. 3. Splitting of a state vector containing burnt and unburnt 
cells into partial vectors with only fuel or ashes. The necessary 
ghost cells at the artificial boundaries are calculated by zeroeth 
order extrapolation. 



n+1 




n+1 




Fig. 4. Determination of the average unburnt interface area 
fraction (3 for two different cases. As can be seen, simply taking 
the average of old and new time level does not always produce 
the correct result. 



and the corner values ■D, r ,i+i/2,j+i/2 are determined by up- 
winding. Depending on the sign of D x at the corner, the ap- 
propriate discrete derivative of G is chosen; if D x is negative, 
one takes (dG/dx) at the right side, and vice versa. Now all 
quantities needed in eq. ( |37| ) are known. 

Because of the discrete nature of the grid, it is in most 
cases impossible to satisfy condition ( f38| ) exactly; therefore a 
re-initialization step is required for the complete implementa- 
tion also. This is done in exactly the same fashion as described 
in section [3.1.2[ 

3.2.4. Flux-Splitting 

In order to compute the total fluxes across a mixed cell inter- 
face, it is necessary to solve the Riemann problems for burnt 
and unburnt states separately. To achieve this, each grid vector 
is splitted into a sequence of completely burnt and unburnt par- 
tial vectors. During this process, artificial boundaries are cre- 
ated at the front locatio n for which boundary con ditions must 
be specified. Following Smiljanovski et al. (1997), this is done 
by zeroeth order extrapolation of the cells at the boundary (see 
Fig. H for illustration). The PPM algorithm implemented in 
PROMETHEUS is then used to calculate the hydrodynamical 
fluxes for the partial vectors. 

Now eq. ( |T(i| ) is applied to compose the total fluxes. How- 
ever, it is in many cases insufficient to use the unburnt interface 
fraction (3 at the beginning of the time step in this formula, es- 
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pecially when the flame enters or leaves a cell during the time 
step. Therefore we calculate the average of (3 over the time step 
(see also Fig. |j): 



7 passive implementation 



1 

At 



(3dt 



(40) 



(41) 



The composed flux then reads 

F = f3F u + (1 - 0)F b . 

3.2.5. Source terms 



The amount of matter consumed by the flame in a mixed cell 
during a time step is given by 



Am = J As u p u dt, (42) 

where A denotes the flame surface in this cell. For the x-sweep 
in a directional splitting scheme one obtains 



n 2 x As u p u dt. 



(43) 



The factor r? x is introduced by the projection of the flame 
on the y-axis (or on the yz-plane in three dimensions) and by 
the projection of the burning speed on the x-axis. In the i-th 
cell, the ratio of the projected flame surface and the surface of 
a cell interface is approximately given by |/3,:+i/2 — 
Thus one obtains for the source terms 



Ashes/, 



At p u . 



A X 

Ax % pi 

AXpuei^ = -AI AsheSii 

Ae to t,i = gAXAshes.i- 

4. Numerical tests 



S u ,i n x,i0i+l/2 ~ A-l/2) 



(44) 

(45) 
(46) 



A set of testcases was run with both of the implementations pre- 
sented above to determine the ability of the numerical schemes 
to represent thermonuclear flames. Our main criteria were the 
reproduction of a given burning velocity and the isotropy of the 
front propagation. Additionally, we investigated the behaviour 
of the algorithms for complex situations, like the merging of 
two flame kernels and cusp formation in a sinusoidally per- 
turbed flame. 

At t = 0, the thermodynamical state of the unburnt matter 
was characterized by p u — 5 • 10 8 g/cm 3 , T u = 5 • 10 8 K, and 
Xizc,u — Xi6 u — 0.5. The energy release for the fusion to 
56 Ni is q = 7 • 10 17 erg/g, and the burning speed s u was set to 
3 • 10 7 cm/s. 

For all tests we used a cartesian grid with a cell size of 
1.5 • 10 6 cm. 




5.0-10 
4.0-10 



7 complete implementation 



0.0 0.2 0.4 0.6 0.8 1.0 

t[s] 



0.0 0.2 0.4 0.6 0.8 1.0 
t[s] 



Fig. 5. Time dependent position of a planar flame propagating 
in positive x-direction for both implementations of the front 
tracking algorithm. The lines indicate the theoretically pre- 
dicted behaviour. 



4.1. ID flames 

In a first test we investigated a planar flame propagating in 
positive x-direction with reflecting boundaries at the left, top 
and bottom edges of the computational domain and an outflow 
boundary to the right. The grid consisted of 128x4 cells. Under 
these circumstances the material behind the front should be at 
rest and the absolute front velocity with respect to the grid is 
expected to be 



Sb := s u p u /pb, 



(47) 



which corresponds to about 4.4 ■ 10 7 cm/s for our initial condi- 
tions. 

As can be seen in Fig. |[ the agreement of simulation 
and predictions is excellent for the complete implementation, 
whereas the passive implementation underestimates the flame 
velocity by about 20%. Figs. || and (7] show the profiles of tem- 
perature, nickel concentration, velocity and density for both al- 
gorithms at t = Is. Again, the complete implementation gives 
exactly the expected results: two constant states that are sepa- 
rated by a mixed cell. With the exception of the fluid velocity, 
the picture is nearly the same for the passive implementation; 
here the transition is smeared out over about three grid cells 
by PPM. The velocity profile shows strong oscillations in this 
case; one also notes that the average fluid motion in the unburnt 
material is noticeably slower than for the complete implemen- 
tation. 

All the deviations in the passive approach can be explained 
by the fact that the flame is not advected with the speed of 
the unburnt matter as postulated in eq. (Q), but by the average 
speed in the burning cells. Depending on whether the flame just 
entered the cell or is about to leave it, this quantity is closer 
to the unburnt resp. the burnt velocity but never reaches the 
desired v u . As a consequence, the flame propagates too slowly 
and at a non-uniform speed, thereby causing fluctuations in the 
velocity field. 

To further investigate this behaviour of the passive imple- 
mentation, two additional tests with p u = 3 • 10 9 g/cm 3 and 
p u = 5 • 10 7 g/cm 3 were performed. In these cases the flame 
propagation speed was underestimated by 14% and 28%, re- 
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Fig. 6. Planar flame propagation test: Temperature, nickel concentration, fluid velocity and density at t —Is for the passive 
implementation. 



1*10 10 
8*10 9 

E 6>1 ° 9 

H 4*10 9 



2*10" 



temperature 



X 



1.5*10 



1.0*10' - 



2*10 7 4*10 7 6*10 7 8*10 7 1*10 8 

x[cm] 

fluid velocity 



5.0*10 



-5.0*10' 




2*10 7 4*10 7 6*10 7 8*10 7 1*10 8 

x[cm] 



1.2 
1.0 
0.8 
0.6 
0.4 

0.2 
0.0 



Ni 







6.0*10° 
5.0*10 8 

4.0*10 8 

S 8 
3.0*10° 

°- 2.0*10 8 

1.0*10 8 




2*10 7 4*10 7 6*10 7 8*10 7 1*10 8 

x[cm] 

density 



2*10 7 4*10 7 6*10 7 8*10 7 1*10 8 

x[cm] 



Fig. 7. Planar flame propagation test: Temperature, nickel concentration, fluid velocity and density at t =ls for the complete 
implementation. 



8 



M. Reinecke et al.: A new model for deflagration fronts in reactive fluids 



spectively. Since the error grows roughly proportionally with 
the density jump, these observations support our interpretation. 
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Fig. 8. Snapshots of the front geometry for circular flame prop- 
agating outwards. The dashed lines represent exact circles and 
have been added to allow easier judgement of the flame geom- 
etry. 
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Fig. 9. Evolution of the front geometry for two merging circular 
flames. The time difference between subsequent snapshots is 
0.1s (from inside to outside). 



However, for the special case of turbulent burning in the in- 
terior of white dwarfs, these seemingly large errors can be tol- 
erated: firstly, the velocity jump across the flame is quite small 
compared to the burning velocity; secondly, our model for the 
turbulent burning speed is based on dimensional analysis and 
therefore s u itself could carry an uncertainty much larger than 
the 28% mentioned above. 

A first order correction for the underestimation of the burn- 
ing speed in the passive implementation can be done in a quite 
straightforward way for this concrete physical problem and will 
be incorporated in future versions of the code. 

4.2. 2D flames 

4.2. 1 . One circular flame 

To test the isotropy of both algorithms, the propagation of an 
initially circular flame was simulated on a grid of 50x50 cells 
with outflow boundaries; some snapshots of the front geometry 
are shown in Fig. ||. While deviations from the circle shape do 
occur, they are sufficiently small for both implementations. The 
difference in the flame propagation speed is still present and 
nearly of the same size as in the one-dimensional simulation. 
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Fig. 10. Evolution of the front geometry for a sinusoidally per- 
turbed flame. The time difference between subsequent snap- 
shots is 0. Is (from left to right). 

4.2.2. Two merging circular flames 

On the same grid as in the simulation above, the merging of 
two circular flame kernels was investigated to demonstrate the 
ability of the level set approach to handle topological changes. 
As the results indicate, the formation of a single front happens 
smoothly and without numerical difficulties (see Fig. ||). The 
slight deformation of the fronts before the merging can be ex- 
plained by the interaction of the velocity fields generated by 
both flames. 

4.2.3. Perturbed planar flame 

Fig. [To] shows the temporal evolution of a sinusoidally per- 
turbed flame propagating in positive x-direction. As expected, 
the trailing part of the front becomes narrower until a cusp is 
formed; afterwards, the flame geometry remains practically un- 
changed. The short vertical section of the flame that can be seen 
in the right panel of Fig. [lO|is an artifact of the rather poor res- 
olution: since the (expected) cusps are located exactly at the 
y-position of the cell centers and the level set is stored at the 
cell corners, they cannot be seen in this discretization. 

4.3. Sensitivity of the reconstruction equations 

The results of all tests described above show that both imple- 
mentations of the level set method can be used to model tur- 
bulent thermonuclear combustion in Type la supernovae. Since 
the complete version is more accurate, it would be the method 
of choice. Unfortunately, however, it has turned out that the 
straightforward implementation of the reconstruction as de- 
scribed above leads to numerical difficulties when applied to 
the "real" situation in a white dwarf including density and pres- 
sure gradients and gravitational forces. In our supernova sim- 
ulations with the complete implementation, after a few time 
steps the reconstruction in many mixed cells failed because the 
internal energy of the unburnt state reached values for which 
the equation of state is undefined: 

&i,u < e E os(pu, T = OK, Xf ue i) (48) 

Discretization errors in the input values are the most likely 
reason for this divergence. To test the reaction of the recon- 
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Fig. 11. Reconstructed temperature of the unburnt material for 
varying deviations in a. Below a/ao < 0.98, the reconstruc- 
tion fails. 



by Miemeyer & Hillebrandt (1995b). Our results show a 
flame which is perturbed due to Rayleigh-Taylor- and Kelvin- 
Helmholtz-instabilities on all scales down to a few grid cells 
(see Fig. |l2|). An extensive discussion of this simulation as well 
as simulations with other initial conditions can be found in Rei- 
necke et al. (1998). 

5.2. Chemical hydrogen combustion 

The complete implementation has already been successfully 
used to model turbulent flame fronts in lean hydrogen-air mix- 
tures. Fig. |l3| shows the merging of three flame parts in a mix- 
ture of 15% hydrogen in air in a box with an outflow boundary 
to the right and reflecting boundaries elsewhere. Since small 
disturbances are amplified by material diffusion in the case of 
hydrogen flames, the burning speed was modified depending 
on the curvature of the front. 



struction algorithm on such uncertainties, the following exper- 
iment was performed: 

From a given pre- and post-front state that exactly fulfill 
the Rankine-Hugoniot jump conditions a mixed state is syn- 
thesized according to eqs. (0)-(^) for an ao equal to 0.5. Then 
a reconstruction is tried for this mixed state, but for a slightly 
different a (i.e. for an a with some uncertainty). In Fig. [l]], 
the reconstructed temperature of the unburnt material is plot- 
ted against the introduced error in a. It can be easily seen that 
for a/ao < 0.98 a reconstruction of pre- and post-front states 
becomes impossible. For highly curved fronts, as they are ex- 
pected in Type la supernovae, the deviations of a from the ex- 
act value can become much higher than that, because a is ob- 
tained for a piecewise linear approximation of the front. At first 
glance, one would expect an improvement if the front geometry 
was modeled with higher accuracy, e.g. by approximation with 
quadrics. But in this case, other problems appear: the number 
of different topological configurations in a cell explodes, and, 
most importantly, there is no way to define the normal n that is 
required by the reconstruction equations. 

Because of these numerical problems, we have not yet been 
able to simulate Type la supernovae with the complete imple- 
mentation of the front tracking scheme. An investigation of the 
properties of the reconstruction equations and, if possible, cre- 
ation of a more robust system is subject of future work. How- 
ever, introducing just an artificial viscosity to limit the curva- 
ture of the flame front may be an easy way to stabilize the nu- 
merical scheme. 



6. Conclusions 

We have presented a numerical model to describe deflagra- 
tion fronts with a reaction zone much thinner than the cells 
of the computational grid. In contrast to the currently favoured 
method for astrophysical simulations ( Khokhlov 1993 ), our ap- 
proach provides a considerably sharper transition from fuel to 
ashes, thereby allowing the growth of hydrodynamical insta- 
bilities on smaller scales and generally the evolution of small 
features in the flame. 

Two different implementations of the model have been de- 
veloped and tested; for simple initial conditions, both versions 
produce results acceptable for our needs. However, because 
of the mentioned numerical problems the complete implemen- 
tation cannot be employed for supernova simulations without 
modification. 

In addition to modeling flames, the level set method de- 
scribed in this paper can also be used for tracking contact dis- 
continuities with only minor modifications; therefore, any ap- 
plication in astrophysical hydrodynamics dealing with one of 
these phenomena might benefit from this numerical scheme. 
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5. Applications 

5.1. Type la supernovae 

The passive implementation of the level set method has been 
used to model the turbulent flame front in the early stages 
of Type la supernova explosions. To allow direct compari- 
son with the reaction-diffusion model, our initial conditions 
were chosen as similar as possible to the simulations done 
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Fig. 12. Temporal evolution of front geometry and velocity field after igniting a single, circular bubble near the center of the 
white dwarf. Note that the scales change from snapshot to snapshot. 




Fig. 13. Chemical combustion in a lean hydrogen-air mixture on a grid consisting of lOOx 140 cells, with reflecting borders at 
the left, top and bottom and a flowout boundary at the right. Due to the high material diffusion of hydrogen, little disturbances 
are amplified. This was modeled by using a curvature-dependent flame speed. 



12 



M. Reinecke et al.: A new model for deflagration fronts in reactive fluids 



References 

Broyden, C. G.: 1965, Mathematics of Computation 19, 577 
Fryxell, B. A., Miiller, E., Arnett, W. D.: 1989, MPA Preprint 
449 

Khokhlov, A. M.: 1993, ApJ 419, L77 

Mulder, W., Osher, S., Sethian, J. A.: 1992, JCP 100, 209 

Niemeyer, J. C: 1994, Turbulente thermonukleare Brennfron- 

ten in WeiBen Zwergen, Master's thesis, Max-Planck-Institut 

fur Astrophysik, Garching 
Niemeyer, J. C: 1995, On the Propagation of Thermonuclear 

Flames in Type la Supernovae, Ph.D. thesis, Max-Planck- 

Institut fur Astrophysik, Garching 
Niemeyer, J. C, Hillebrandt, W.: 1995a, ApJ 452, 769 
Niemeyer, J. C, Hillebrandt, W.: 1995b, ApJ 452, 779 
Nomoto, K., Thielemann, F. K., Yokoi, K.: 1984, ApJ 286, 644 
Osher, S., Sethian, J. A.: 1988, JCP 79, 12 
Press, W. H., Teukolsky, S. A., Vetterling, W. T, Flannery, 

B. P.: 1992, Numerical Recipes in C, Cambridge: Cambridge 

University Press 
Reinecke, M. A., Hillebrandt, W., Niemeyer, J. C: 1998, MPA 

Preprint 1 122b, accepted by A&A 
Sethian, J. A.: 1996, Level Set Methods, Cambridge: Cam- 
bridge University Press 
Smiljanovski, V., Moser, V., Klein, R.: 1997, Comb. Theory 

and Modeling 1, 183 
Sussman, M., Smereka, P., Osher, S.: 1994, JCP 114, 146 
Timmes, F. X., Woosley, S. E.: 1992, ApJ 396, 649 



